Velocities of transmission eigenchannels and diffusion

The diffusion model is used to calculate both the time-averaged flow of particles in stochastic media and the propagation of waves averaged over ensembles of disordered static configurations. For classical waves exciting static disordered samples, such as a layer of paint or a tissue sample, the flux transmitted through the sample may be dramatically enhanced or suppressed relative to predictions of diffusion theory when the sample is excited by a waveform corresponding to a transmission eigenchannel. Even so, it is widely assumed that the velocity of waves is irretrievably randomized in scattering media. Here we demonstrate in microwave measurements and numerical simulations that the statistics of velocity of different transmission eigenchannels are distinct and remains so on all length scales and are identical on the incident and output surfaces. The interplay between eigenchannel velocities and transmission eigenvalues determines the energy density within the medium, the diffusion coefficient, and the dynamics of propagation. The diffusion coefficient and all scattering parameters, including the scattering mean free path, oscillate with the width of the sample as the number and shape of the propagating channels in the medium change.

is symmetric about the center of the sample.Further, since the deviation from linearity of the energy density is small, as seen I (b), we include only the leading order, quadratic, correction to the constant gradient of intensity within the sample.This gives an excellent fit.The parameters of the fit for this sample with  = 1.078 are given in Supplementry Note 11.(b) The energy density within the sample is seen to fall nearly linearly with depth for samples with length ~1.

Supplementary Notes 1. Recursive Green's function simulations
Simulations of electromagnetic wave propagation are carried out on model rectangular samples shown schematically in Supplementary Fig. 2. Samples of length  and width  are comprised of a lattice of squares with sides  = /2 and random values of the dielectric constant, , drawn from a rectangular distribution [1 − Δ, 1 + Δ] with Δ = 0.3.The sample is open on the left and right and bounded on top and bottom by perfect reflectors.Thus, field excited from the open boundaries of the sample may be expressed in terms of the waveguide modes of the empty waveguide.
We employ the recursive Green's function method [3][4][5] to find the field in the k th column for a source on the left hand side of the sample.The sample is divided into K columns of width , labeled 1, 2, …, k, … K. To find the field in each column of the sample, we first consider the field in an isolated column with random dielectric constant in each element within an otherwise empty waveguide.The field in the isolated k th column is found by solving the homogeneous wave equation (  −   )Ψ = 0, where Ψ is the electric field,   = ∇ 2 +  2 , and   is the eigenvalue of the wave equation.We refer to   as a Hamiltonian because the wave equation is discretized in analogy to the tight-binding Hamiltonian used in electronic systems 3 .
We seek to solve the equation �  −   (,  ′ )�  (,  ′ ) = δ( −  ′ ) for a point source.We are interested in the causal solution for the Green's function,   = lim →0 + (  +  −   ) −1 .The total Hamiltonian of the system is: Here,   is the Hamiltonian of the left lead,   is the Hamiltonian of the right lead,   is the Hamiltonian of the k th column, and   is the interaction term between the neighboring k th and l th columns.We begin with the surface Green's function of the left lead and proceed from left to right, connecting each column to the subsystem to its left.The coupling between a k th column and the subsystem to its left is expressed via the Dyson equation,  =  0 +  0 . is the Green's functions of the connected columns from 1 to k.This is the matrix of solutions of the Hamiltonian submatrix from 1 to k.  0 is the matrix of Green's functions of the connected columns from 1 to k-1 and the isolated column k, for a subsystem with the Hamiltonian and  is the interaction terms between the two, given by This gives the recursion equation where the superscript  indicates the connected subsystem starting at the left side,   is the Green's function of the connected subsystem at column k beginning from column l, and A similar recursion relation is found by iterating from right to left.The left-to-right and right-to-left parts are then combined to obtain the result,  1 = [1 −     ,+1  +1,+1   +1, ] −1  1  , which allows us to calculate the Green's functions for all columns from 1 to K. The Green's functions are used to calculate the energy density, TM, transmission eigenvalues and EVs, and the energy density in each TE.

Scaling of eigenchannel velocities for 𝑵𝑵 = 𝟔𝟔𝟔𝟔
In contrast to the exponential decay of transmission eigenvalues with length in samples with  = 8, the EVs asymptotically approach a constant value with increasing length, as seen in Fig. 2(a,b).Most EVs approach an asymptotic value at lengths shorter than the localization length of approximately ℓ s , but the EV of the transmission eigenchannels with the lowest transmission approaches its asymptotic value at a length which is comparable to ℓ s .To check whether the approach of EVs to their asymptotic values is related to the onset of multiple scattering or of localization, we carried out simulations in samples with the same disorder of ∆ = 0.3 but with  = 64.The localization length is then much longer than the scattering mean free path.The results in Supplementary Fig. 3 show that the EV of low transmission eigenchannels close approach their asymptotic values for lengths decidedly shorter than ℓ s .This indicates that the evolution of values of EVs with length is related to ℓ s and not ℓ s .Thus, the EV do not change in this sample in the crossover to Anderson localization and their distinct values are not related to Anderson localization.

Comparison of microwave and optical measurements of transmission
The nature of transmission and energy density within a random medium depends upon the excitation.When the TM can be measured, it is possible to excite the medium with the sum of all incident waveguide modes, as is the case in the present paper.This occurs naturally in measurements of electrical conductance.it is possible to excite Lower values of EVs for channels with smaller   is consistent with measurements of a drop in transmission of a laser beam directed at a random dielectric slab with increasing angle of incidence, which corresponds to a reduced normal component of the velocity of the light [26].The transmission through the , where  p is the distance traveled by the beam before the direction light is randomized, and  ′ is the angle of refraction within the sample.
The derivative of energy density is constant within a random medium and  b is clearly defined when the flux is the same in all incident channels.More generally, as for the case of a laser beam incident on a random slab, the energy density will be a weighted superposition of the eigenchannels and will produce a different spatial profile of energy density, even deep within a random medium.

Relationship between incident and reflected transmission eigenchannels
For a dissipationless quasi-1D medium with reciprocity, the scattering matrix can be expressed as  = �   ′   ′ �, where  and  are the reflection and transmission matrices for a wave incident from the left and the corresponding matrices for excitation from the right are primed.
Reciprocity implies   = , while the conservation of flux in a lossless medium implies  †  = 1.The polar decomposition of the scattering matrix [ ,] gives where   and   are unitary matrices.Here, With these relations, and the knowledge of ,   and   , we can obtain  ′ , ,  ′ and, as thus the full scattering matrix.In the main text, which deals primarily with the TM, we have indicated   by  and   by .Thus, for the transmission matrix, the incident singular vector is   , and the transmitted singular vector is   .For the reflection matrix, the incident eigenvector is   , the reflected singular vector is −  * .
The velocity of the n th transmitted eigenchannel is v , = ∑   � ,, �

Energy density at the input surface
Here we show that the average energy density at the input surface is equal to the sum of averages of the energy densities of the incident and reflected waves,   (0) =  ,i (0) +  ,r (0).This is a consequence of the proportionality of the reflected transmission eigenchannel and the complex conjugate of the incident TE, which was shown in Appendix D. Since the demonstration is based on relationships in single configurations, we distinguish between a property of a single configuration and the average over configurations by using a bracket for the latter.

Normalization of energy density profile by its spatial average
Energy density profiles for ensembles with different values of  are shown in Fig. 6(a).The energy density profile () is normalized by its average over configuration and position.This average is obtained by noting that the average energy densities excited in each transmission eigenchannel on the opposite side of the sample for excitation from the left,   (), and the right, . Thus, the sum of energy density at the input for excitation in all channels on both sides of the sample is . Using Eq. ( 2), gives, (0 . Together with the result above that  , =  v + , this gives (0) =  , /2 6,7 .Since the spatial average of the energy density excited by all channels from either the left or right are equal,

Expression for nonlocal diffusion coefficient
Substituting the expression for  in Eq. ( 6),  = , into the expression for  in Eq. (4b) gives Straightforward manipulation gives, which is the result of Eq. ( 7).

Determining the scattering mean free path
Since the scale of disorder in the scattering model considered in simulations is much smaller than the wavelength,  = /2, scattering is nearly isotropic, and resonances do not affect the energy velocity.We therefore expect the scattering mean free time,  s , to be nearly independent of the direction of propagation of the wave.This is confirmed in the analysis below, which shows that the scattering mean free times for different incident waveguide modes, which have different transverse velocities, are identical.The scattering mean free path ℓ  in a sample with Δ = 0.3,  = 64 and width  = 64.5  0 2 is found from the decay of the coherent flux vs. ballistic time delay of the waveguide modes.The scaling of the coherent flux of waveguide modes, |〈  ()〉| 2 , is shown in Supplementary Fig. 4(a).|〈  ()〉| 2 is the ensemble average of the flux of the  th waveguide mode at the output of a random sample of length .The waveguide modes with higher index, , which have smaller group velocities, v w , decay more rapidly.The decay of |〈  〉| 2 vs. coherent time delay /v w , is shown in Supplementary Fig. 4(b).The decay times of the coherent flux for the m th waveguide mode is 8 ,  s, = ℓ , v  .Since all the plots collapse to a single curve and fall exponentially to give a single mean free scattering time  s , and the speed of the wave in the medium is c, the scattering mean free path is ℓ  =  s = 16.0 m.

Pulse propagation
Under steady-state illumination, (, ) can be expressed via Fick's first law in terms of the parameters of the TEs,   and v  , as in Eq. (4a).In the time domain, (, ) is given by Fick's second law as the ratio of the partial derivative of concentration with time and the Laplacian of the concentration.In samples thicker than five times the mean free path, the transmitted intensity in a diffusive medium falls exponentially in time once higher order diffusion modes have decayed, with a decay rate 1/ D =  2 /( + 2 b ) 2 9-11,10 .This is the decay rate of the quasinormal modes or resonances of the medium and equals the linewidth of quasi-normal modes 12,13 .
The transmission due to an incident Gaussian pulse is obtained by multiplying field transmission spectra with the Gaussian spectrum of the envelope for the field pulse and Fourier transforming into the time domain.The ensemble average of the transmitted intensity for the sample for which steady-state results are shown in Fig. 4(a) with  = 64 and  = 7ℓ  , at which   0 � reaches its peak value of 0.91, is shown in Fig. 4(a).The transmitted pulse decays exponential with decay time  D = 0.68 μs ± 0.01.Beyond a delay of 5 , the decay rate slows due to the increasing relative weight of longer lived quasi-normal modes [12][13][14] .The diffusion coefficient given by  = Though the diffusion coefficients found in steady-state and in time-domain simulations agree, still  falls 10% below  0 � .Several factors contribute to this difference.At the length at which () reaches its peak value, the diffusion coefficient is already lowered below the bare diffusion coefficient due to weak localization 12,13,15 .The bare diffusion coefficient might be obtained by extrapolating the linear region of ()/ 0 � in Fig. 4(a) back to  = 0.This gives, (0)/ 0 � = 0.95.In the diffusive regime, the length dependent diffusion coefficient, (), falls linearly with length from the bare diffusion coefficient with the fractional reduction in () proportional to 1  ℎ ~ℎ / T because of the increasing probability of a path looping back upon a coherence length of the path with increasing sample thickness [10].Similarly, the time dependent diffusion coefficient proportional to the decay rate of ln〈()〉 falls linearly with delay time, with the drop in decay rate proportional to / T 12,13 .
Another reason that / 0 � may fall below unity is that  0 � = 1 2 ℓ  may be larger than the bare diffusion coefficient  0 = 1 2 v E ℓ 16 .There are two countervailing factors, ℓ may be larger than ℓ  , while v E may be smaller than .These corrections are likely to be small since  =  0 /2, so that the individual scattering elements are too small to support resonances so that v E ~, and scattering might be nearly isotropic with ℓ~ℓ  .However, it should be noted that the prevailing assumption that scattering is isotropic and is independent of sample dimensions, which leads to the classical result for the diffusion coefficient, is shown here not to hold for wave diffusion in bounded media.

Scaling of diffusion coefficient and its factors
The scale-dependent diffusion coefficient (, ) relative to  0 � = 1 2 ℓ  and the factors of the diffusion coefficient  b (, ) and v T (, ), as given in Eq. (7), are given in Supplementary Fig. 6. Results for samples with disorder ∆ = 0.3 and for  = 8 and  = 64 are shown, respectively, in the top and bottom rows over the full range of  for each  and for a wide range of /ℓ  .The variation with width is smaller for larger , but is still substantial for  = 64.

Breakdown of diffusion in energy density profile
We find that energy density within the sample falls linearly with depth for sample with lengths nearly equal to the localization length  =  and  = 1.Here we consider the variation of energy density in an ensemble of samples with  = 1.078.The sample width is  = ( +

7 . 7 |
Departure from linearity of () near the localization threshold Supplementary Fig. Departure from linearity of () near the localization threshold.(a) The magnitude of the gradient of the energy density is greatest at the center of the sample.Since the position-dependent diffusion coefficient and the derivative of the energy density depend upon the distance from the nearest boundary, (()/) v +

(𝐿𝐿+2𝑧𝑧 b ) 2 𝜋𝜋 2 1 2
D11 , which is equal to 2.16 × 10 9 ± 0.02 m 2 /s, obtained with  b taken as its average value over the spectrum for the incident Gaussian pulse,  b = 11.4 .The average scattering mean free path over the spectrum of the incident pulse of ℓ  = 16.0 ± 0.15 m, gives  0 � = ℓ  = 2.40 × 10 9 m 2 /s and / 0 � = 0.90 ± 0.015, in agreement with the state-diffusion coefficient for this sample, as seen in Fig.4(a).The uncertainty in / 0 � is due to different values of the slope of the curve in the inset obtained for different time ranges.